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ABSTRACT 


A theoretical engineering thermal analysis of man-made “permafrost 
islands" as offshore drill sites is conducted by developing a two- 
dimensional numerical model from the heat conduction equation. The 
model is formulated by solving an equivalent = principle 
based on the finite element method. Phase transformation due to 
reeciting and thawing is assumed to occur isothermally to promote 
Solution economy. Solutions in the time domain are obtained by 
the Crank-Nicolson method. 

Boundary conditions may be arbitrarily imposed or simulated at the 
air-surface interface. Surface boundary conditions are simulated by 
using a one-dimensional surface energy exchange model which is coupled 
to the two-dimensional model. Local meteorological data is employed 
for the surface energy simulation. Successful tests were conducted 
that compared the numerical model solutions with exact solutions to 
freezing front propagation. Also, comparisons were made between model 
results and in situ measurements of soil temperatures at Fairbanks, 
Alaska. Finally, the model was applied to the problem of offshore 


“permafrost island" drill sites. 
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Chapter 1 


INTRODUCTION 


1.1 Alternatives to Conventional Offshore Drilling Structures 


Recent studies indicate that offshore regions in the Beaufort and 
Chukchi Seas may contain extensive oil and gas reserves. Some geologists 
Beit eve that offshore areas tend to be more productive than adjacent 
land areas. If this is true for arctic waters, the implications are 
Staggering considering the recent petroleum discoveries along the arctic 
coast. It can be expected that exploratory offshore drilling in the 
Beaufort Sea will commence sometime within the forseeable future with 
initial probes restricted close to shore. 

Offshore drilling experience in hostile arctic environments has 
been limited to Cook Inlet, Alaska. Although annual growth of sea ice 
in Cook Inlet is only 3 feet, design wind, wave, and earthquake forces 
on offshore structures have been found to be negligible when compared 
with expected lateral ice forces. Obviously, the coastal environment 
in the Beaufort Sea is much more hazardous to ocean structures since 
annual sea ice reaches thicknesses of 7 to 9 feet and is subject to 
continuous motion due to the driving forces of wind and current. 

Other dangers include the possibility of pressure ridges, ice islands, 
or multi year ice striking a structure. 

Assuming an adequate knowledge of critical design loading due to 
ice, the application of known engineering techniques for offshore drill- 
ing structures could probably be physically implemented but, at an 


economically prohibitive expense. The arctic thermal regime may permit 


consideration of two additional alternatives for development of offshore 


drill sites: 

1. Construction of man-made permafrost islands from dredged sea- 
floor materials that are frozen into an icy matrix strong 
enough to resist lateral forces exerted by annual sea ice 
movement. 


2. Construction of structurally reinforced man-made ice islands 


grounded at the desired location. 


1.2 Scope and Objective 


This report is addressed to the first alternative. The primary 
problem, thermal analysis of nearshore man-made permafrost islands, is 
considered. The objective is to develop a versatile engineering tool 
_based on known thermodynamic theory and make appropriate assumptions to 
facilitate engineering application. A mathematical model is formulated 
for computer solution utilizing the finite element method. Numerical 
solutions are compared with analytical solutions prior to application 
to example permafrost island problems. This report is restricted to a 
theoretical study relying on available data for soil thermal properties 
and available historical records for meteorological conditions. The 
results of several long-term simulations are presented. They are 
considered to be of value to engineers desiring design information on 
equilibrium freezing front propagation in various soil structures 


subjected to offshore arctic meteorological conditions. 





Chapter 2 


FINITE ELEMENT FORMULATION OF TRANSIENT HEAT CONDUCTION WITH PHASE CHANGE 


2.1 Freezing and Thawing of Soils 

Heat transfer with phase change due to freezing or thawing is the 
principal thermal process of interest in applied arctic soil engineering 
problems. When attempting to mathematically quantify or simulate the 
physics of these processes, the engineer must assess what level of model 
sophistication is necessary in the context of the proposed application. 
The primary factors affecting his decision are the extent to which the 
physical properties of the soil system are known and the Aieey of 
information available for determining external boundary conditions. 

As a general rule, the physical parameters of the soil system will 
dictate how the phase change problem should be solved. Three alterna- 
tive methods are considered for modeling the latent heat of fusion 
problem for saturated soils: 

A. Phenomenonological 

‘B. Thermodynamic 

C. Analytical - isothermal 


Each method and its inherent limitations is briefly discussed. 


2.1.1 Phenomenonological Modeling 

It is well known (Keune and Hoekstra, 1967) that certain moist or 
Saturated soils do not undergo isothermal phase change. Lukianov and 
Golovko (1957) proposed a phenomenonological equation to describe such 


a process in a soil-water-ice system: 





93> _d¢ 
Ox Kx) = Foe Q) 


where x = distance from the air-soil interface (ft) 
t = time (day) 
@ = temperature (°F) 
K = K (x,t,¢) = thermal conductivity (BTU/ft-day-°F) 

F = apparent volumetric heat capacity (BTU/ft>-°F) = C - L dG/d¢ 


C, = C, («,t,$) = volumetric heat capacity of the soil-water-ice 
system 


G(?) = Ice content (lb/ft?) 


L = latent heat of fusion (BTU/1b) 


Equation (1) mathematically describes the freezing-thawing phenome- 
hon in a one-dimensional zone of finite thickness. The freezing of a 
small differential element is accompanied by a temperature dependent 
release of latent heat which decreases the time rate of change of 
internal energy if net heat extraction from the element remains invariant. 
The phenomenonological approach is equivalent to introducing an "apparent 
volumetric heat capacity" as defined above. 

Although phenomenonological modeling of phase change remains 
analytically intractable, it is possible to numerically solve equation (1). 
Nakano and Brown (1972) developed a one-dimensional finite difference 
scheme from equation (1) and satisfactorily simulated heat conduction 
in tundra soils near Barrow, Alaska. It is important to stress that 
the success of any phenomenonological model is heavily dependent on 


the following factors: 





1. Unfrozen soil water content as a function of temperature must 
be .known or experimentally determined. 

2. Numerical simulation of phase change requires small time incre- 
ments to avoid solution instabilities caused by the injection 


or removal of latent heats. 


3. Applied boundary conditions must be known on a micro-scale. 


2.1.2 Thermodynamic Modeling 

Thermodynamic modeling of soil phase change also incorporates 
equation (1) to predict soil column temperatures as a function of time 
but does not rely on experimentally determined soil freezing curves to 
Simulate the temperature dependent evolution of latent heat. Recent 
research (Koopmans and Miller, 1966) has shown that many soil moisture 
characteristic curves and soil freezing curves are similar. Since a 
sOil moisture characteristic curve defines the relationship between 
water content and soil pore water pressure, detailed calorimetric 
analyses to determine unfrozen water content in soils as a function of 
temperature can be avoided if temperature-pore water pressure relation- 
ships during phase change are known. 

It is possible to calculate soil water pore pressure in soils under- 
going phase transition by assuming that pure ice and water can only 
co-exist in thermodynamic equilibrium as defined by the Clausius - 


Clapeyron equation for a bulk fluid (Sears, 1964) 


LdT 
dp= Tay) (2) 
where L = latent heat of fusion 
T = existing absolute temperature 





v" = specific volume of ice 
v' = specific volume of water 
dp = incremental change in pressure 


It is sometimes assumed that the Clausius - Clapeyron equation applies 
to unsaturated soils and that ice pressures in the soil pores are at 
atmospheric pressure. It is probable that this equation does not apply 
to saturated soils without modification. 

A typical numerical procedure for thermodynamic simulation of soil 
freezing 1s to compute a soil temperature profile from equation (1), use 
this to obtain soil water pore pressure by numerically integrating 
equation (2), and enter the soil moisture characteristics curve to 
determine the change in water content. The net change in water content 
_ over a given time increment determines the magnitude of latent heat 
that must be released in the following time step. A new or updated 
apparent volumetric heat capacity is computed to continue the temporal 
solution of equation (1). This method of discretely modeling phase 
transformation requires the utilization of small time increments to 
avoid unstable oscillations in computed temperatures. Furthermore, 


soil moisture characteristic curves must be experimentally determined. 


2.1.3 Analytical Modeling of Isothermal Phase Change 


Analytical solutions commonly employed in the engineering analysis 
of phase change problems are simplified versions of the theory of ice 
growth introduced by Franz Neumann in the 1860's. Examples are Stefan's 
(1889) and Berggren's (1943) applications of the Neumann theory to ice 


growth and frost penetration problems. Since any theoretical formulation 





of phase change (commonly referred to as the "problem of Stefan") is 
nonlinear, classical solution techniques are not applicable (Carslaw 
and Jaegar, 1959). 

Because of the mathematical complexities involwed, “the praéticrng 
engineer is usually forced to seek a published analytical (exact) 
solution that represents the best approximation to his physical problem. 
A significant disadvantage of mathematically exact solutions is the 
limitation imposed by simplifying assumptions and solution space geo- 
metry. To analytically solve the isothermal phase change problem, it 
is usually necessary to assume a semi-infinite isotropic medium initially 
at a uniform temperature (i.e., simplified boundary conditions are 
chosen so that the governing equation of state can be solved). The 
imposed boundary condition is a sudden change in temperature that subse- 
quently induces phase change in the medium. It is the restrictive 
imposition of simple boundary conditions that severely limits the 
general application of exact solutions. 

Several Russian investigators have been successful in analytically 
modeling more realistic boundary conditions. Kolesnikov (1946) presented 
a derivation which incorporated the major hydrologic energy exchanges 
affecting the growth of sea ice, and more recently, Portnov (1962) 
developed an exact solution for isothermal phase change with time 
dependent boundary temperatures. 

Because of the relative paucity of exact solutions and the complexi- 
ty of modern engineering problems, most research efforts have been 
directed toward achieving approximate methods of solution. The most 


effective have been finite difference and finite element methods which 





rely on analytical solutions as a means of testing numerical techniques 
prior to simulating physical problems that remain analytically insoluble. 
Hence the primary importance of known theoretical solutions is that of 
investigating the accuracy of a numerical procedure whose validity has 


not been theoretically established. 


2.2 Derivation of a Two-Dimensional Finite Element Model 

An approximate numerical solution to the parabolic differential 
equation of transient heat conduction in two dimensions is developed by 
employing the finite element technique based on the so-called Rayleigh- 
Ritz formulation. Phase change is assumed to occur isothermally to 
avoid previously discussed difficulties common to thermodynamic and 


phenomenonological methods of solving the latent heat problem. 


2.2.1 Variational Finite Element Formulation 

The variational functional approach for the finite element formula- 
tion of transient heat conduction relies principally upon the theoretical 
contributions of Biot (1970) who introduced a unified variational 
analysis of heat transfer analogous to the energy theorems of classical 
mechanics. Given the partial differential equation of two-dimensional 


transient heat conduction in Cartesian coordinates 


Bry 8b) 4 Ory 2 = eee 3 
5x Ka) * oy Kay? + G = pea (3) 
where © = temperature 
K, and are the thermal conductivities in the x and y 


direttions, respectively 
G = G(x,y,t), a variable heat source or sink 


Oc = volumetric heat capacity 





t = time 
the transition to an equivalent variational statement is sought by 
applying an extension of the fundamental lemma of variational calculus 
(Weinstock, 1952). The theorem is applied by noting that if in a given 


time interval, the functional 


oan ad 96 
7 EGY Ose mo ramen (4) 


is to be minimized over a bounded region R, a necessary and sufficient 
condition is that at any given instant of time the unknown function 


o(x,y,t) satisfy the Euler equation 


a of d of Or 
3x F(ab7exy’ * By BCBe7Byy? ~ 39 ~ © - 
subject to ¢ obeying the same boundary conditions of equation (3). The 
problem is to obtain a functional which when substituted into equation 
(5S) reduces to equation (3). 
Assuming so remains invariant, an equivalent variational formula- 
tion of equation (3) is obtained if the functional 


dg 


7 She Ux (Be en? + KyG we) + 2(pca- - G)o}dxdy (6) 


can be minimized for all admissible states of $. For equation (6) to 

be a true minimizing principle, the second variation of the functional 
with respect to the unknown state variable must be greater than zero. It 
can be shown that this is indeed true for transient heat conduction 
(Myers, 1971). It is less restrictive to require the functional to be 
Stationary (i.e., only the first variation is zero) and demonstrate 


that the numerical solution will converge to an exact solution. 
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Equation (6) is formally solved by discretizing the domain R into 
numerous subregions or "elements" and assuming that the temperature 


Ch 


distribution within each element can be approximated by an n~ order 


polynomial. In this report, a linear polynomial is selected to 
approximate the state variable in a two-dimensional finite element mesh 
composed of arbitrarily shaped triangular elements. 

Upon substitution of a linear temperature distribution shape 
function into equation (6), the first variation of the functional is 
evaluated by applying Leibnitz's rule for differentiation of integrals. 
Thermal properties are assumed to be constant within each element although 
abrupt changes in properties between elements are permitted. Summation 


of the integrated contributions of all elements produces a system of 


Simultaneous linear first order differential equations of the form 


[K]{o} + [c] {<2} = {G} (7) 
where [K] = system conduction matrix 
[C] = system capacitance matrix 
{G} = system load column matrix 
{o} = nodal temperature column matrix 


The system conduction and capacitance matrices are banded and symmetric, 


Significantly reducing computer memory requirements for storing arrays. 


2.2.2 Incorporation of Isothermal Phase Change 


Simulation of heat conduction and phase transformation in two 
dimensions involves a minor modification of the variational formulation 
of transient heat conduction. As previously discussed, the solution 


domain is discretized by a finite element mesh composed of triangular 
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_ elements. The mesh geometry is stored in computer memory by reading in 
nodal points, their respective global coordinates, thermal properties 
of each element (frozen and unfrozen), and moisture content of each 
element. After this has been accomplished, numerical simulation of 
phase change may commence. = 

Consider a common node point "A" in a segment of a two-dimensional 
finite element mesh (Figure 1) whose temperature is dependent upon any 
arbitrary boundary condition imposed along S. Phase change is not 
permitted to occur at node "A" until a requisite amount of latent heat 
has been exhausted or added. The prevailing quantity of latent heat 
for node "A'' is determined by adding one third the weight of material 
(subject to phase transformation) of each contributing element to a 
nodal latent heat accumulator array assigned to node ''A", 

As an example, consider the cooling and freezing of node "A". 
After each time increment that a solution is being generated, the temp- 
erature of node "A'' is scanned to determine whether or not it has 
dropped below a prescribed point where phase transformation is assumed 
to occur (TFAZ). If the node has not been previously frozen and the 
computed temperature has dropped below TFAZ, then the quantity of 
latent heat evolved during the previously computed time step becomes: 

DELQ = C,, (TFAZ - Computed Temperature) * Volume (8) 
where C,, = a weighted volumetric heat capacity based on the heat 

capacities of the soil-water mixture assigned to node "A" 
that must undergo phase change. 
Volume = volume of soil-water assigned to node ''A" 
The quantity DELQ is subtracted from the latent heat accumulator for 


node '"'A'' and if more latent heat must still be exhausted, the computed 





Fig. | 


Typical 


NODE "A" 





Finite Element Mesh 
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temperature at the node is brought back to TFAZ prior to obtaining a 
new solution at the end of the next time increment. This is analogous 
to simulating phase change as an isothermal process. The time rate of 
latent heat generation is governed solely by the solution of the heat 
conduction equation and remains insensitive to large time steps. 

If the required amount or more than the required amount of latent 
heat has been extracted, node "'A'' and the volume of soil-water assigned 
to it are considered to be frozen. The thermal properties of the 
elements common to node "'A'' are then updated based on the volumetric 
proportions of soil-water and soil-ice present in the respective 
elements. If an excess amount of latent heat is removed, the residual 
is used to calculate how far below TFAZ the temperature at node "A" 
should be. A weighted C¢ is used for this computation. All weighted 
volumetric heat capacities (Cy and Cr) are computed from the thermal 
properties of elements common to nodal points undergoing phase change. 

The thawing process is handled in a similar manner. Provisions are 
made to monitor which nodes are frozen or thawed so as to avoid re- 
freezing frozen nodes. Boundary nodes remain unaffected by the checking 
routines that scan nodal temperatures. After phase transformation at 
a node occurs, the latent heat accumulator assigned to that node is 
recomputed. Recomputation permits the simulation of long-term cyclic 


freezing and thawing. 


2.2.3 Imposition of Boundary Conditions 
The variational functional expressed by equation (6) automatically 


incorporates boundary conditions of either specified temperature or 
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zero heat flux. Inclusion of additional modes of heat transfer that 
may occur over various portions of the boundary surface, shown in 
Figure (1), 1s accomplished by adding appropriate surface integrals 

to equation (6) which upon minimization yield the desired boundary 
conditions. The surface integrals for specified heat flux, convection, 


and radiation heat transfer across the boundary are 


Be 2B. 


S atds + 1s K(o? - 2bap)ds + of (Legs? - eqdhods (9) 
B B. > 

where q = specified heat flux across the boundary 

d = boundary node temperature 

d.. = ambient. temperature near the boundary 

K = convective heat transfer coefficient 

Oo = Stefan-Boltzmann constant 

Ep = Emissivity of the boundary surface 

Ex = Emissivity of the ambient surroundings 
Outgoing heat fluxes are considered positive. 

The combination of equations (6) and (9) and their subsequent 
minimization produce a generalized equivalent variational statement 
describing transient heat conduction with its associated boundary 
conditions. Substantiating mathematical details are provided in 
Appendix A. Minimization of surface boundary integrals produces terms 
that are combined with the system conduction and load matrices of 
equation (7). Nonlinear radiation heat transfer terms give rise to 
solution difficulties since unknown temperatures are introduced into the 


System matrices. Methods of circumventing this problem include lineari- 
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zation or modeling net radiation exchange as a specified heat flux. 


These methods are discussed in greater detail in section 3.2.2. 


2.2.4 Solution in the Time Domain 
Linearization or removal of nonlinear terms from the system 
matrices allows preservation of the previously discussed system of 


ordinary differential equations represented by 
[K]{o} + [c] {2%} = {a} (10) 
Ot 


and greatly simplifies the time domain solution. After a considerable 
amount of experimentation with various numerical time domain solution 
techniques, the Crank-Nicolson method was selected since it provided 
the most satisfactory results in terms of optimizing time step sizes, 
of providing acceptable solution accuracy, and of minimizing numerically 
induced oscillations. 

The Crank-Nicolson method utilizes the arithmetic mean of the 
derivative of the state variable at the beginning and end of each time 


step to move the solution ahead in time. Thus 


(g}V*2 = (oY + SECtghY + (94) (11) 
where At = time step size 
v = a reference point in time 
Upon premultiplication by the capacitance matrix and substitution of 


equation (10), equation (11) becomes 
({K] + Sch) {e}’** = GEfc) - [k]){o}” + 2{6) (12) 


Since {o}" is the initial or previously computed temperature 


distribution and the system matrices [K], [C], and {G} remain constant 
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in each time step, equation (12) can be further simplified to the 
matrix equation 

(w]{o}"*2 = {ZH (13) 
where [W] = a symmetric matrix in band form. 
Instead of employing a time consuming matrix inversion algorithm, the 
solutions, {¢}¥t!, are computed by the Gaussian elimination method which 
exploits the symmetric and banded structure of the [W] matrix. 
2.3 Comparison of Finite Element Solutions with Analytical Solutions 

Involving Phase Change 

Because the method formulated to solve the latent heat problem is 
a numerical approximation, comparison of finite element solutions with 
several exact solutions to simple isothermal phase change problems is 
justifiable. Two problems are considered: one in which thermal para- 
“meters remain unaffected by phase change and secondly, one where phase 
transformation produces a significant variation between frozen and 
unfrozen thermal properties. 

.2.3.1 Simulation of Isothermal Phase Change with No Change in 
Thermal Parameters 

As a first example, the problem of the solidification of a homo- 
geneous slab of finite thickness is considered. An exact solution 
discussed by Allen and Severn (1952) is used for comparison purposes. 
In Figure 2, a two-dimensional finite element mesh for analyzing the 
problem is shown. Figures 3 and 4 summarize the finite element results 
obtained from two different time domain solutions: implicit and Crank- 
Nicolson. Good agreement between both numerical solutions and the exact 


solution is evident. Note that the Crank-Nicolson time domain solution, 
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Figure 2. 


Figures 3 and 4. 


Solidification of a slab of liquid: finite element 
mesh. 


Distance between nodes l. and 19 = L = 60. AX = L/6 = 10 
K = 1; pe = /Y2 (thermal properties remain constant) 
Uniform initial temperature = 45° 

Temperature of phase transformation = 45° 

Latent heat of transformation = 70.26 pc/unit volume 


The entire perimeter with the exception X = L is 
thermally insulated. 


The end X = L@t ¢=0 is dropped to and subSe- 
quently maintained at a temperature of 0°. 
Solidification of a slab of liquid: progress of 


freezing front through the slab. 


Inset diagram illustrates nomenclature in each 
rigure. 


The freezing front is denoted by line F at a typical 
time t, and the frozen zone is shown shaded. 


Exact solution for the equation of the solidification 
locws 1S: 


2_ kK 
(X-L)* = —2(1.07)t 


as obtained from Jones' theory (Allen and Severn, 
1952) 
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Figure 4, permits similar solution accuracy with time steps twice as 
eve as those employed in the implicit time domain solution, Figure 3. 

This example demonstrates the accuracy of finite element discretiza- 
tion of the freezing process. Judicious selection of time steps and 
finer mesh subdivision will improve solution accuracy with increased 
computer simulation time. A much more rigorous and realistic test is 
to consider the following theoretical problem where substantial changes 
in thermal properties of the conducting medium occur upon phase change. 

2.3.2 Simulation of the Neumann Theory of Freezing Front Propaga- 

tion in a Saturated Soil Column 

The application of the Neumann theory to the prediction of frost 
penetration depth in soils is well known. A transcendental equation is 
solved to obtain a constant of proportionality whose magnitude is 
dependent upon the thermal properties (frozen and unfrozen) of a semi- 
infinite soil column, initial uniform temperature distribution, and 
imposed initial constant boundary temperature. This constant yields the 
exact solution for freezing front propagation when multiplied by the 
square root of elapsed time. 

Comparison of a one-dimensional finite element simulation with a 
Neumann solution for freezing front propagation is shown in Figure 5S. 
The Neumann solution was extracted from an example cited by Jumikis 
(1966). The freezing front locus is plotted versus the square root of 
time to demonstrate the stability of the finite element solution and 
enable computation of the numerical constant of proportionality. The 
theoretical constant of proportionality is .575 ft/day 1/2 and the finite 


element solution constant of proportionality is .56l ft-day 2/2, With 
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Figure 5. Comparison of a variational and analytical solution 
to the Neumann problem. 


A. 


Physical constants: saturated sand at a void 
ratio of .50; specific gravity = 2.65. 


Thermal parameters: 
1. Unfrozen: 


Ka 


C, = 42.7 BTU/ft®°F 


25.7 BTU/ft-°F-day 


ul 


2. Frozen: 


Ke 
Cf 


32.2 BIU/ft-°F-day 


29.3 BTU/ft2°F 


Initial and imposed boundary conditions: 


Initial uniform temperature distribution of 
36°F. A constant boundary temperature of 14°F 
is maintained. Isothermal phase transformation 
is assumed to occur at 32°F. 


Solutions: 


Analytical and variational solutions are: 


.575 vt (theoretical) 


x 


.561 Yt (finite element) 


x 


where x = frost penetration depth (feet) 


t = time (days) 

Variational solution was obtained by employing 
a one-dimensional finite element mesh. Spacing 
between nodes was one foot with numerical 
computations conducted in .5 day increments. 
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such a small discrepancy in constants of proportionality, it would 

take an elapsed time of 15 years wEberE the error between computed and 
theoretical frost penetration depth exceeded one foot. At that time the 
theoretical freezing front would have penetrated 42.6 feet below the 


ground surface. 
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Chapter 3 


SIMULATING AIR-GROUND INTERFACE ENERGY EXCHANGES 


3.1 Methods of Approximating Heat Transfer at the Air-Ground Interface 

The thermal processes occurring at the interface between a soil 

System and the atmosphere are often too complex for an exact descrip- 
tion; yet, their very existence is crucial to conducting an analysis 

of the thermal regime. The dominating processes are absorbed incoming 
solar radiation, net long wave radiation exchange between the atmosphere 
and surface boundary, turbulent heat transfer, evaporation, snow and 
vegetative cover effects, and the soil temperature gradient at the air- 
ground interface. There are three alternative methods that may be 
employed to simulate the net effect of these processes: 

1. Impose a soil surface boundary temperature based on known 
ambient air temperature or in situ temperature measurements 
of the soil surface. 

2. Compute an equivalent temperature acting across a known 
thermal resistance to produce a net flux into the soil surface. 
With this method soil surface temperatures are computed rather 
than imposed. 

3. Mathematically model each of the individual thermal processes 
as precisely as possible and solve the resulting system of 
nonlinear equations simultaneously with equation (7). 

The method formulated in this report for simulating soil surface 

energy exchanges is a hybrid combination of alternatives 2 and 3. 


Individual thermal processes are modeled on a macroscopic scale, con- 
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verted to an equivalent temperature, and coupled to a variable turbulent 
heat transfer coefficient to simulate heat fluxes transmitted across the 
air-surface interface. This approach generates surface temperatures 


and eliminates the necessity of solving a system of nonlinear equations. 


3.2 Modeling the ieee aes peeeaees afernecar Exchange 
From a macroscopic viewpoint an energy balance at the air-ground 
interface can be approximated by the expression 
Qsoil + Qswr + Qniwr + Sconv - Qevap = 9 (14) 
where Q.54) = heat flux into the soil column 
Qcwr = Short wave radiation component 
Qniwr = het long wave radiation exchange term 
Qcony = turbulent heat transfer component 
Qevap = evaporative heat loss term 
and where all heat fluxes are expressed in BIU/ft2-day. 
Prior to discussing how equation (14) is incorporated into the variation- 
al analog for transient heat conduction, the methods employed to deter- 


mine each of the above energy components are briefly discussed. 


3.2.1 Short Wave Radiation 

Incident short wave radiation is determined by 

1. Applying CRREL nomographs (Gerdel et.al., 1954) for computing 

| available solar energy as a function of time and latitude. 

2. Modifying available solar energy with empirical relationships 
(Eagleson, 1970) which take into account altitude of cloud 


base and per cent cloud cover. 
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3. Selecting appropriate surface albedos from Wechsler and Glaser 


(1966). 


3.2.2 Long Wave Radiation Exchange 

As discussed in section 2.2.3, the net radiant heat emission term 
is nonlinear. If the surface emissivities of two radiating bodies are 
the same, a linearization technique discussed by Mettler (1966) and 
Pivovarov (1973) could be applied to provide 

Quiwe 7 40&¢3,(6-$,.) (15) 
where all terms are as defined in section 2.2.3. 

Since the emissivity of a soil or snow surface is different than 
that of the atmosphere, it is easier to compute net long wave radiation 
exchange from the Stefan-Boltzmann law and apply net radiation heat 
transfer as a discrete flux across the boundary surface. An "apparent 
atmospheric emissivity" is assigned from empirical formulas (Eagleson, 
1970) which incorporate per cent cloud cover, cloud base altitude, and 
atmospheric water vapor pressure as influencing variables. 

Inaccuracies are introduced because previously computed boundary 
surface temperatures are used to determine net radiative heat transfer 
in the next time increment of the transient heat conduction problem. 
This is preferable to formulating an exact variational statement because 


computer solution by iteration is avoided. 


Seno Evaporative Heat Losses 


Evaporative heat transfer is not directly simulated but is estimated 
from data and applied as an outgoing flux across the surface boundary. 


The magnitude of evaporative flux for a given region is determined by 
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analyzing the potential evapotranspiration (PET) or actual evapotranspira- 
tion (AET) losses from data given by Patric and Black (1968). A time 
dependent evaporative heat flux function 1S constructed with a temporal 
shape closely approximating that of available solar energy. Its inte- 
grated area equals the selected PET or AET. In this report, it is 
desired to simulate the effect of evaporation from a saturated soil 
system; hence, PET is employed to compute an evaporative heat flux 
function. Evaporation during the winter months is assumed to be neglig- 


ible. 


3.2.4 Turbulent Heat Transfer 


Turbulent heat transfer is modeled according to the relationship 


= K(¢ - 9 (16) 


Qe onv surface a 


a turbulent heat transfer coefficient that is a 
function of wind velocity, surface roughness, and 
atmospheric stability 


where K 


Oeameace = surface temperature 


daiy = ambient air temperature 
The theoretical and experimental contributions of Scott (1964, with 
subsequent verification in 1969) are used to compute a turbulent heat 


transfer coefficient that 1S allowed to vary at the end of each time 


increment or when meteorological parameters change. 


3e2e5 wSnow Covem Effects 

It is not desired to simulate the hydrologic morphology of a snow 
cover over the soil surface but to observe the gross insulating effect 
on a seasonal basis. In this respect heat conduction through an average 


winter snow cover is modeled using the previously discussed energy 
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exchanges. Snow cover density is estimated from data presented by 
Biello (1969) and thermal conductivity from curves discussed by Mellor 
(1964). 

Minor sophistications are included to permit degradation of snow 
surface albedo when the ambient air Penne snes rises above 32°F. 
Simulation of heat conduction through the snow cover is terminated 
when the computed snow surface temperature is greater than or equal to 


32°F. 


3.3 Finite Element Simulation of Surface Energy Exchange 

Having determined individual heat flux components, net energy 
exchange across a two-dimensional surface is simulated by coupling 
one-dimensional elements having a thermal conductivity equal to the 
turbulent heat transfer coefficient and having a volumetric heat 
capacity of zero to the two-dimensional elements representing the 


boundary surface. An equivalent temperature is calculated from 


— Qewr + Qniwr a Qevan 17 

oe _ i K (17) 
where Peq = equivalent temperature 
domb = ambient air temperature 


K = turbulent heat transfer coefficient 
and other terms remain as previously defined. The computed temperature 
from equation (17) is then imposed at the end of each one-dimensional 
element. ' 
The imposition of equivalent temperatures at discrete nodes produces 


a net uniform flux (satisfying equation 14) across the two-dimensional 


boundary surface once they have been properly inserted into the system 
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matrix equations. Although substantiating mathematical details are 
omitted for the sake of brevity, note that if the equivalent tempera- 
ture is equal to the ambient air temperature, only turbulent heat 
transfer is simulated. Figure (6) is a schematic presentation of the 


concept. 


3.4 Simulation of Surface Energy Exchange under Fairbanks Field Conditions 
Numerical simulations of surface energy exchange were conducted 

for Fairbanks, Alaska summer and winter climatic conditions to permit 

evaluation of several solution techniques. Numerical instabilities 

encountered in modeling winter snow cover effects led to the development 

of the equivalent temperature approach for handling thermal processes 


at the air-surface interface. 


3.4.1 Surface Energy Simulation: 1-10 July 1972 

Using local meteorological data for Fairbanks, summer surface 
temperatures of exposed silt were simulated over a ten day period. Silt 
thermal properties were computed from the experimental data of Kerstens 
(1949). A clear day outgoing evaporative heat flux was determined from 
Fairbanks evapotranspiration loss data (Patric and Black, 1968) and 
further modified by a factor equal to the ratio of actual solar energy 
received to theoretical available clear day solar energy. Input data 
and results are shown in Table 1. Computed ground surface temperatures 


were approximately 8°F warmer than the ambient air temperature. 
t 


3.4.2 Surface Energy Simulation: 20-29 January 1973 
A ten day winter simulation of heat conduction in a soil column 


with an 18 inch snow cover was conducted to allow comparison of 





ou 


NET FLUX INTO SOIL 


pie dd Edd 





SYMBOLS. 
Peq = Equivalent Temperature 
ow-O0 = One-Dimensional Surface 
Energy Exchange Element 
K = Turbulent Heat Transfer Coefficient 
pe = Volumetric Heat Capacity=O 
Fig. 6 


Modeling Surface Energy Exchanges 





32 
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measured and computed temperature profiles at the end of the 10 day 
period. Of secondary interest were the computed snow surface tempera- 
tures. Average Fairbanks January snow density and thermal properties 
were determined from values presented by Wheeler (1973). In this 
Simulation, solar and long wave radiation fluxes were placed directly 
across the snow surface without employing an equivalent temperature 

and turbulent heat transfer coefficient. Solution instabilities caused 
by day to day fluctuations in meteorological parameters were encountered 
when all surface energy balance terms were simulated in such a manner. 
Available options to solve the stability problem were those of utiliz- 
ing very small simulation time increments or applying the equivalent 
temperature concept. The latter was chosen. 

Results of the winter simulation are shown in Table 2. Ten 
elements were used to simulate the snow cover and six elements to model 
a 45 inch soil column. A lower soil boundary temperature of 32°F was 
imposed since field observations indicated a mean January position of 
the groundwater table at a depth of 45 inches below the ground surface. 


Calculated and measured results were in good agreement. 
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Table 2 
Winter Simulation (20- 29 January 1973) 
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Chapter 4 


TWO-DIMENSIONAL PERMAFROST ISLAND THERMAL ANALYSIS 


4.1 Statement of the Problem 
A theoretical engineering analysis was conducted of the thermal 
feasibility of man-made permafrost islands. Two different problems 
were considered: 
1. The thermal degradation of an island drill site built on 
location during the winter months and completely frozen 
prior to spring breakup. 
2. The propagation of a winter freeze front through an island 
constructed from dredged seafloor sediments during the summer. 
Simulated island geometry was a right cylinder with a vertical 
height of 25 feet and radius of 24 feet immersed in 20 feet of water. 
Computer storage limitations prohibited investigation of larger island 
geometries. Two-dimensional finite element discretization of a vertical 
cross section of the island is shown in Figure (7). Symmetry in the 
transient heat conduction solution was assumed, thus, only the left 
half of the island is shown. 
Daily meteorological data from Barter Island, Alaska, was utilized 
to simulate surface boundary temperatures. Barter Island is located 
2 miles from the north shore of the Alaskan mainland in the Beaufort 
Sea. The climate is determined by the surrounding ocean environment 
which prevents extremely low winter temperatures and limits the warming 


effects of continuous daylight during the summer months. 
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4.2 First Thermal Analysis: Prefrozen Drill Site 

The first thermal analysis conducted assumed a prefrozen drill 
site built during the winter at an initial uniform temperature equal to 
the ambient air temperature (15°F). Thermal simulation was for the 
period 1 May 1970 through 31 August 1971. ~ 

The upper 4 feet of the island consisted of saturated gravel with 
a porosity of .45 and the remaining 21 feet were saturated silt with a 
porosity of .5. Soil thermal properties were determined from Kersten's 
(1949) report. Gravel in the upper part of the island that thawed 
during the summer was assumed to drain to a moisture content of 10 per 
cent of gravel unit dry weight. Draining would deter thawing, and 
gravel thermal properties were modified accordingly. Subsequent winter 
freezing was assumed to occur isothermally at 100 per cent soil satura- 


tion at 32°F. 


4.2.1 Applied Boundary Conditions 


To conserve use of computer time, upper surface boundary conditions 
were estimated by conducting a one-dimensional finite element surface 
energy exchange simulation. A vertical cross section of the island was 
employed. Barter Island meteorological data was provided in 24 hour 
increments, and daily surface temperatures were computed for use in the 
two-dimensional model. Outgoing summer evaporative heat flux was 
estimated from Figure (8). The presence of a 15 inch mean winter snow 
cover at a density of .33 g/cm was simulated for the period 1 October 
1970 to 15 May 1971. Monthly averages of simulated energy terms are 


shown in Figure (9), where average net flux across the boundary surface 
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| is derived as a residual term. In Figure (10), simulated average 
ronthly surface temperatures and snow-ground interface temperatures 
are shown along with average measured ambient air temperatures for 
comparison purposes. 

The effect of a .02°F per foot geothermal gradient acting across 
a sea bottom of frozen silty clay was simulated and found to be neglig- 
ible. For this reason, the bottom of the two-dimensional island was 
considered to be nonconducting. 

A mean annual seawater temperature of 28.5°F was applied along 


the completely immersed vertical sides of the island. 


4.2.2 Simulation Results 

Two-dimensional computer simulation results are summarized in 
Figures (11) through (14). These figures show computed isotherms for 
the 1970-71 winter and following summer for the cases of mean seasonal 
and negligible winter snow cover. The second case assumes that snow is 
removed to promote maximum soil cooling. In both simulations, approxi- 
mately 2.5 feet of the surface gravel thawed by the end of the summers 
of 1970 and 1971. To determine porosity effects on thaw, a one- 
dimensional simulation predicted a 5 foot depth of thaw if gravel 


porosity was .30 instead of .50. 


4.3 ere Thermal Analysis: Dredged, Unfrozen Drill Site 

A second two-dimensional thermal analysis was undertaken to 
determine winter freezing front propagation in a dredged artificial 
island with the same geometry and soil structure as the previously 


discussed case. Natural cxtraction of heat was enhanced by assuming 
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that no significant amount of snow would be allowed to accumulate on 
ihe surface of the island. 

Salinity of the seawater in the saturated silt and gravel-was 
assumed to be 15 parts per thousand which depressed the freezing point 
to 30.4°F. Oceanographic measurements (Hufford, 1974) indicate that 
salinities may very along the arctic coast from 3 to 20 parts per 


thousand in the upper 50 feet of water. 


4.3.1 Imposed Boundary Conditions 

As before, a one-dimensional thermal analysis of surface energy 
exchange was accomplished utilizing a saturated silt soil column to 
determine if enough silt would freeze over one winter to justify a two- 
dimensional thermal simulation. Results are shown in Figure (15) where 
depth of freezing front penetration is plotted against degree days of 
frost. It was judged that enough silt did freeze to consider a two- 
dimensional analysis worthwhile. 

It was also recognized that the imposition of a mean annual sea 
temperature of 28.5°F would not be completely representative of imposed 
winter boundary conditions along the vertical sides of the island. For 
this reason, the growth of sea ice was simulated and the temperature 
distribution in the ice was imposed along the upper submerged portion 
of the island considered to be in contact with the ice. Geothermal heat 


flux effects were again assumed to be negligible. 


4.3.2 Simulation Results 
Two-dimensional simulation results are shown in Figure (16). 


Approximately the upper 13 feet of the island froze during the first 
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winter. To assist in predicting the elapsed time for the entire 

‘s sland to freeze, a one-dimensional simulation incorporating geothermal 
heat flux was run with the results shown in Figure (17). It was 
estimated that the island would be completely frozen after 5.5 years. 
Extrapolation of results indicated an equilibrium frost penetration 
depth of 27 feet would be reached in a saturated silt soil structure 
with a 4 foot saturated gravel over-burden. 

A second one-dimensional simulation using saturated gravel pre- 
dicted that a gravel island of the same dimensions would be completely 
frozen in 3.5 years. Computed equilibrium freeze front penetration in 
gravel subjected to Barter Island climatological conditions was 48 feet 
after an elapsed time of 18 years. Simulation results are shown in 


Figure (18). 


4.4 Discussion of Results 

Both thermal simulations demonstrated that surface degradation of 
island drill sites due to summer thawing would be less than 3 feet, 
which is quite minimal. The primary explanation for this is that 
average Barter Island summer cloud cover is approximately eight-tenths, 
and consequently the magnitude of transmitted short wave radiation is 
low. As an example, average simulated transmitted solar heat flux, 
from Figure (9), during the month of June was only about 50 per cent 
of the available theoretical clear day solar energy. 

The convergence of computed isotherms (Figures (11) and (13)) near 


the air-water interface indicated that a large quantity of heat was 
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extracted from the surrounding water, transmitted through the island, 
and exhausted to the atmosphere .during the winter. 

Of considerable concern is the intrusion of heat along the vertical 
face of the island during the summer months as shown by the position of 
the 29°F isotherm at the end of the 1971 summer (Figure (14)). The 
danger of thawing can be eliminated by constructing prefrozen island 
drill sites with low salinity water. Dredged islands would not present 
a major problem since brine cells should migrate towards the freezing 
front during the winter and be rejected into the unfrozen silt. The 


high salinity silt would then become diluted during the summer months. 
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Chapter 5 


CONCLUSIONS AND RECOMMENDATIONS 


5.1 Conclusions 

The objective of this report was to formulate and develop a viable 
numerical model designed to solve arctic engineering problems concerning 
transient heat conduction accompanied by isothermal phase change. The 
first phase of the work involved comparison of the numerical model with 
several analytical solutions to phase transformation problems. Numeri- 
cal results and analytical solutions were found to be in good agreement. 

A second aspect of this research dealt with the simulation of 
surface boundary temperatures by employing meteorological parameters 
that are readily available for most geographic locations of interest. 
It must be emphasized that the theoretical and experimental results of 
numerous investigators were relied upon to evaluate the individual 
hydrologic energy exchange terms discussed in Chapter 3. Simulated 
energy terms and surface boundary temperatures were the results of 
simplified approximations to complex thermal processes. The numerical 
simulations of surface energy exchange under Fairbanks summer and winter 
climatic conditions assisted in assessing the simplifying assumptions 
employed. Of particular interest were the small discrepancies between 
Simulated and measured soil column temperatures at the end of the 10 
day Fairbanks winter simulation, considering that ambient air tempera- 
ture fluctuated from -39°F to 16°F. 

The final phase of the work applied the thermal model to an 


engineering analysis of silt and gravel offshore drill sites. Based 
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on the thermal simulation results, the following conclusions are made 
subject to Barter Island, Alaska, climatological conditions: 
1. Under natural freezing conditions with negligible snow cover, 
equilibrium frost penetration depth for saturated silt with 
a 4 foot gravel overburden is approximately 27 feet and for 
Saturated gravel, 48 feet. 

2. Summer thaw of saturated frozen gravel at a porosity of .45 
can be expected to be 3 feet and at a porosity of .30 can 
be expected to be 5 feet. 

3. Dredged island drill Sites are thermally feasible in shallow 

water depths. 

Expanding on the above conclusions, the depth of water a dredged 
island could be built in would depend on a number of variables such as 
the magnitude of lateral sea ice forces and the manner in which the 
island is dredged. A possible method of constructing a dredged island 
is to place the silt or gravel in an expendable metal torus perforated 
to relieve excess pore water pressures as the island freezes. Use of 
a torus would prevent loss of material due to wave erosion and hasten 
freezing in contrast to dredged soils deposited at a small natural 
angle of repose. Since both silt and gravel freeze faster than seawater, 
the metal torus would not be destroyed unles the compressive strength 
of the frozen materials was exceeded or the frictional resistance of 
the entire structure was less than the maximum expected lateral force 
due to sea ice movement. 

Lateral resistance to sea ice could be enhanced by driving piles 


at the center of the torus and by using sea bed foundation design methods 
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currently employed on gravity structures in the North Sea. The 
 eeeiion of an inclined rubble breakwater around the island would 
permit rafting of sea ice and reduce the horizontal component of ice 
forces. Furthermore, where drill sites are built in extremely shallow 
water, it becomes possible for seafloor sediments beneath the sites to 


freeze and increase resistance to lateral movement. 


5.2 Recommendations for Future Research 

Thermal analysis of permafrost island drill sites is no longer zi 
problem of academic interest alone. Several drill sites have been 
constructed and occupied by Imperial 0il Limited of Canada in very 
shallow portions of the Mackenzie Delta (Riley, 1974), and plans are 
being developed to construct a dredged island drill site in 25 feet of 
water during the summer of 1974. As arctic offshore exploration 
progresses into deeper water, much larger ice forces will be encountered 
relative to the nearshore areas where ice mobility is somewhat limited. 
It will soon become necessary to develop methods of constructing ina 
relatively short time span massive structures capable of withstanding 
substantial ice sheet pressures. An economically feasible structure 
worthy of design consideration would be a grounded ice or ice-aggregate 
island. 

Before any large scale prototype ice island is built, a considerable 
amount of engineering research effort will be directed toward solving 
the thermal problem of ice accretion. An optimum water thickness will 
be sought which will be easy to apply, will minimize thermal degradation 


of ice already produced, and when repeatedly applied, will yield the 
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desired thickness of ice within the available time constraints. Of 
secondary interest will be the problem of determining the volume of 
seafloor sediments capable of freezing beneath the base of an ice 
island. This is of value in computing resistance to horizontal forces 
and estimating the size of the island. 

It is recommended that the finite element method be considered as 
a research tool for these problems. With minor modifications to the 
numerical model formulated in this report, useful design solutions to 
the ice accretion problem could be generated. No modifications are 
necessary to solve the seafloor sediment freezing problem. In addition 
to the thermal problems, it is also recommended that laboratory tests 
be conducted to evaluate freezing front boundary shear strengths of 


Beaufort Sea sediments at probable future offshore drill sites. 
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Appendix A 


MATHEMATICAL DETAILS 








Variational Solution to the Heat Conduction Equation 


It is stated in Chapter 2 that, upon minimization, the generalized 


functional 


x = SILK, GE)? + KGS? + 20058 - G)ohdxay + Lads 


(A.1) 
+ Lh(G" - 2hab)ds + of end? - egdlo)ds 
where cis that portion of the boundary surface where values of $ are 
not imposed. 

satisfies the heat conduction equation and the associated boundary 
conditions of radiative, convective, and specified heat flux. The 
transition from equation (A.1) to the heat conduction equation will be 
demonstrated. 

Assume at a given instant in time ae invariant and let the 
first term in equation (A.1) be represented by any arbitrary functional 
of the form 

LhE(XY 5b, 0x5 by )dxdy (A.2) 


Taking an arbitrary variation of f and its derivatives, produces 


= IALGGSO) + Ge80,) + Geddy) + fadods 


aby 
(A.3) 
4 
+ £h($5o - $.56)ds + of (e_o" 5h - Exnbodh)ds 
Because the 6 and <_ operators commute and a minimum of the first 
variation is desired, the previous equation becomes 
of of 3 of 9 
Ox = Lplagse + 3d Fx 659) + 36,3 py of exeY 
(A.4) 


+ LIq8d + h(45d - 4,80) + o(€,978h - €n0269)]ds = 0 








A-2 


The last two terms in the first integral of equation (A.4) may be 


integrated by applying Green's theorem of integral calculus such that 


SESE 5(60) + By Sb) daxdy = oleae ) + xs 57) Jaxdy 
(A.5) 
- £8 af 
é: > Cxs5 + 36, ——)ds 


where 1x and ly are the direction cosines in the x and y directions 
respectively 


Upon substitution equation (A.5) into equation (A.4) we have 


_ of 0 of 0 ,of 
dx = Lpsolse - 3x 36,0 = By Bg, 14xey (A.6) 


f of 4 4 _ 
x + ‘30, a q = h(¢ ca boo) a o(Ep a ExoPoo) Jds = 0 








+ £89 [1x3 


Since 69 may be any arbitrary variation in R and ¢ is not imposed on c, 


application of the fundamental lemma of variational calculus requires 


that in R 

36 - TxGG) - Typ = ° (A.7) 
and along c 

xg + L3G + a+ AC - Ga) + fend - eWhS) = 0 (A.8) 


By setting f equal to the first term of equation (A.1), it can be 
seen that we have an equivalent solution to the heat conduction equation 


if é6x is indeed a minimum. 
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COMPUTER PROGRAM 
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User Instructions 
a e e e . e 
This finite element program solves one and two-dimensional 


transient heat conduction problems involving isothermal phase trans- 
formation. Boundary temperatures may be specified or internally genera- 
ted at the air-surface interface by providing specified meteorological 
parameters. One hundred and fifty elements and 100 nodes are permitted 
With the restriction that the numerical difference between any two 
numbered nodes common to an element not exceed 19; i.e., a maximum 
system bandwidth of 20 is permitted. All nodes must be consecutively 
numbered beginning with one and all elements must remain within the 
first quadrant of the Cartesian coordinate system. If one-dimensional 
elements are used, they are required to be vertically arranged. 

Energy simulation across a two-dimensional surface is accomplished 
by the use of one-dimensional elements as discussed in Chapter 3. The 
one-dimensional surface nodes and elements must be consecutively 
numbered from left to right beginning with one. The remaining two- 
dimensional elements and their respective nodes may be numbered in any 
desirable sequence. Do not specify imposed boundary temperatures on any 
of the one-dimensional nodes used for surface energy simulation. 
Equivalent temperatures are computed and imposed within the program 
itself. During data input assign a zero volumetric heat capacity and 
any ple eiaiane conductivity to the one-dimensional elements. New 
conductivities are computed after meteorological data has been read in. 

Thermal units of BTU/FT/DAY and °F must be employed if the surface 
energy simulation mode of the program is selected. Any system of 
consistent units may be used for imposed boundary value problems with 


the exception that time be expressed in days. 





Data Input 
FIRST CONTROL CARD (613, 5SF10.2) 


A. 


Cc. 1-3 


CC. 4-6 
CC. 7-9 


CC. 10-12 


CC. 13-15 


CC. 16-18 
CC. 19-28 
CC. 29-38 


CC. 39-48 


CC. 49-58 


CC. 59-68 


1 for one-dimensional finite element simulation; 
2 for two-dimensional finite element simulation; 


3 for coupled one and two-dimensional elements for 


surface energy simulation. 

Number of elements. 

Number of specified temperature boundary nodes. 
Number of nodes not held at a specified temperature 
less number of nodes used for surface energy simulation. 
Total number of nodes. 

Number of one-dimensional elements used for surface 
energy simulation. 

Latent heat of transformation per unit weight of 
substance that is subject to phase transformation. 
Temperature at which isothermal phase transformation 
occurs. 

Time step size in DAYS or decimal fractions thereof. 
Desired elapsed time between nodal temperature print- 
outs, expressed in days. 


Total simulation duration in days. 


SECOND CONTROL CARD (15) 


CC. 1-5 


Number of elapsed time segments before new boundary 


or meteorological data is to be read in. 
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THIRD CONTROL CARD (1F4.1); (Delete if surface energy simulation 
not desired) 


CC, 1-4 Average snow depth present. If no snow cover.is being 
Simulated, enter a zero. 

FOURTH CONTROL CARD (1F5.2); (Delete if surface energy simulation 

not desired) 

CC. 1-5 Enter 1 if it is desired to reform system matrices in 
every time increment due to large fluctuations in 


meteorological parameters; otherwise enter a zero. 


FIRST DATA CARD (8F7.2) 

CC. 1-7 List initial temperature of nodes, excluding specified 
temperature boundary nodes, in node number sequence. 

CC. 8-56 Repeat above procedure. Use additional cards as 


necessary. 


SECOND DATA CARD (2F10.2) 
Cc. 1-10 X coordinate of node. 
Cc. 11-20 Y coordinate of node. 


CARDS MUST BE IN NODE NUMBER SEQUENCE! 


THIRD DATA CARD (313, 8F5.1) 


Cc. 1-3 First node number of element. 
CC. 4-6 Second node number of element. 
CC. 7-9 Third node number of element. Enter a zero if the 


element is one-dimensional. 
Node numbers of two-dimensional elements must be read 


in a counter clockwise direction around the element. 








CC, 
cc. 
CC. 
CC. 
cc. 
cc. 


cc. 


cc. 


10-14 
15-19 
20-24 
25-29 
30-34 
35-39 


40-44 


45-49 
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Existing thermal conductivity of element. 
Thawed thermal conductivity of element. 
Frozen thermal conductivity of element. 
Existing volumetric heat capacity of element. 
Frozen volumetric heat capacity of element. 
Thawed volumetric heat capacity of element. 
Per cent moisture content of element. (Per cent of 


element unit weight). 


Element unit weight. 


ALL CARDS MUST BE IN ELEMENT NUMBER SEQUENCE! 


FOURTH DATA CARD (7(13,F7.2)) 


CC. 


cc. 


cc. 


1-3 


4-10 


11-70 


Node number. 
Specified boundary temperature. 


Repeat above sequence. 


DATA MUST BE PROVIDED IN NODE SEQUENCE! 


FIFTH DATA CARD (16,2X,7F8.2); 
not desired) 


CC. 


cc. 


cc. 


CC. 


CC. 


CC. 


1-6 


41-48 


(Delete if surface energy simulation 


Date (day, month, year), i.e., 071271. 
Ambient air temperature (°F). 
Wind velocity in miles per hour. 
Cloud cover in tenths. 

Cloud base altitude in 1000 feet. Enter 34 if 


ceiling is unlimited. 


Relative humidity in per cent. 
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CC. 49-56 Clear day solar heat flux (BTU/ft?-day). 


CC. 57-64 Outgoing evaporative heat flux (BTU/ft*-day). 


J. REPEAT H and I based on the number of times data is to be updated 


or boundary conditions changed. 


Computer Program Listing 


Complete computer program listing and example output are provided 
at the end of this appendix. Printed energy balance terms of evaporative 
heat flux, net long wave radiation exchange, and convective heat flux 
are considered as heat losses from the air-surface interface when 


positive. 





~ 7000 
7025 


7050 FORMAT(1L0X-! 


— me - 


7075 
7100 
7159 
7200 

__ 7250 
7275 
7300 
7350 


~~ 7400 


7450 
_ 7509 
7510 
7520 
7525 
7530 
7550 
7600 


Se ee ee ee 


8000 
8050 


8100 
~ 8150 


8175 


~ 8290 


_ 8225 


8240 
8250 


i 


8275 
8300 


8325 


~~ 8350 
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INTEGER DATE 


FORMAT (61375F10.2) 
FORMAT (15) 


— a a em ee me me ee eee 


DATA ERROR- UNSPECIFIED NODES+ BOUNDARY NODES DO 

1 NOT EQUAL TOTAL NUMBER GF NODES® ) -~ _ ” a 7 oa 
FORMAT (1F4.1) 
FORMAT (8F7.2) | 
FORMAT(2F10.2) © 
FORMAT(313+8F5.1) 
FORMAT(T(E3,F7.2)) 
FORMAT(1642X%,7F8.2) 
FORMAT (LOX,38HDATA FERROR-SYSTEM BAND WIOTH TOO 
FORMAT(LOX,22HNEGATIVE AREA ELEMENT,I5) 

FORMAT( 10X»4HNODE» 2X91592XeL3HFROZEN ON DAY s2XsF722) 
FORMAT (10Xs4HNODE, 2% 515 e2XeL3HTHAWED ON DAY22X,F 702) 
FORMAT(10X,*COMPUTED NODAL TEMPERATURES ON DAY" »2XeF7.2//) 
FORMAT{(7(2Xe"*NODE* ¢3X_* TEMP.*)) a 
FORMAT(2X%e7( *—--—---—------')) 
FORMAT (//¢2X9T( teen maa) 
FORMAT (2X»7( °-------——------ *)s//) 
FORMAT(T(2X%e1493XeF6-2) ) 

FORMAT (10X», "SEMULATION OVER SPECIFIED TINE SPAN HAS BEEN 
1 COMPLETED’ ) ~ 2 ell US eee 
FORMAT (LOXe3Xe’SHVP® 93X09 3Xo *QSOL® 9 3X9 2X "EVELX® 9 3X0 2X_*QLRAD’, 
13X03X»*QCONV? » 3X) 

FORMAT(//410X%-"CGMPUTED SATURATION VAPOR PRESSURE AND ENERGY 
IBALANCE TERMS CN DAY*,FT7.2) 

FORMAT (10X,5F10.2) 

FORMAT(1LOX~e*SNOW COVER HAS MELTED SNOW COVER SIMULATION TERMINATE | 
LD ON DAY® sF7e2e%e'e/elLOXe*COMPUTED NODAL TEMPERATURES ARE * ) 
FORMAT(1LOXe’SIMULATION CALENCAR DATE IS(YRMODY) %,2XeISe//) 
FORMAT (1LOX,*’SOLUTION HAS BECCME UNSTABLE ENERGY SIMULATION HAS 
LBEEN TERMINATED? ) 

FORMAT (LOX,s*COMPUTED EQUIVALENT TEMPERATURE® ,1F6.2¢13X, "AMBIENT 
LAIR TEMPERATURE’ ,1F5.2) ' ae 
FORMAT(//,10%e"THE FOLLOWING NODES ARE PARTIALLY FROZEN *) 
FORMAT(10X,*NODE*,5X»"PERCENT OF ASSIGNED VOLUME FROZEN*,«5X,*ASSIG 
NED VOLUME*) © ' | 2 
FORMAT(LOXs14¢19XeLFL0.2931X» LF 542) 

FORMATCLOX»s*DATA ERROR NO BGUNDARY NODE SPECIFIED. AT LEAST ONE _ 
LBOUNDARY CGNOITICN MUST BE READ IN AS DATA®’) 5 ; 
FORMAT(//,19X%s"NO NUDES UNDERGOING PHASE TRANSITION AT THIS TIME’? 

1//) 

FORMAT(LFS.2) 


—— eee 


a 


a a ee eS 


LARGE) 


Ee et 


me ee ee 


— a mee ee ee er eee — — © 


tee ee 


Po 5 I As a ee i Ge i om 


COMMON 
COMMON 


~ COMMON 


COMMON 
COMMON 


~ COMMON 


CCMMON 
COMMON 


~ COMMON 


COMMON 


/JBLKL/ 
/BLK2/ 
/BLK3/ 
/BLK4/ 
/BLK5/ 
/BLKG6/ 
/BLKT/ 
/BLKB/ 
/BLKIO/ 


R( 100) 


S( 100,20) 


ST( 190,20) 
NOD(15003) | 


X{ 100) 
¥( 100) 
TK (150) 


ROwWC (150) 
/BLKLO/ WC 1LC0,20) 





sr re ee ee eee ee 


oC mee we ee ee 





4... = COMMON /BLKLI/ Z(1CO) | ae _ 5 he, Ve 
COMMON /BLK12Z2/ NBC(59) 


CUMMON /BLKL3/ BCISO) 
DIMENSION RWCF(100),2WCT( LOO) »NKN( LOO) » TOL 100) 2G 100) 
DIMENSICN TXT(1352) 9 TKFL150),RCWCT( 159) sROWCF(1L50) »CMNOD(150),. 
LFLUIO( 150) sFNOO( 150) ePChTRI(150) PUNWAT(150) 
___ OIMENSICN VGL(100) ,FZHET(100) »NFREZ( 100) »WATE( 100) 
DIMENSICN A(3-3) »PEL(3,3),XLCC(3),YLOC(3),SEL(3,3) 
READ(1L,»7000) KODE s NEMS eL BC eo NK eo MM yp NSURF pHEAT y TFAZsOT pDAY END 
READ(1,7025) NTSEG 
=" GGUNN=EAY/O8 9 “Se aass..  «@ 


—=— a 





ee EE ee ee ee ee 





COUNTI=0. 
TIME = 0. 
es > TEAY= 0. °°" ee cic mee ee = a 
ITSEG = 0. 
_ SNOW = 0. 
oS OPT = O. : a a er EEE 


IF «{NSURF) 200,200,100 
100 READ(1,»7075) SNOW 
a a READ( 11,8350) OPT 
200 CONTINUE 
IF (MM~(NK4+L6C) «EQ. 0) GO TO 500 
——— WRITE(3,7C50) 
GO TQ 33900 
500 CONTINUE 


ae ee eee -- SS RS A | 
( ee 


C—----READ INITAL NOOAL TEMPERATURES 
C----— 
READ(L,71C0) "© TO (1) y T= 15MK VO 
00 1000 I=1l,yNK 
| VOL{I) = 0. 
ae AE Ile Oe 
RWCF(I) = 0. 
RWCT(I) = 0. 
TFI(TOCIY .GT. TEAZY GO TO ~~ 800 
NFREZ(I) =1 
GO TO 1C00 
===" —"" 800 NFREZ(I) = 2 
1000 CCNTINUE 


see eee ee ee eee - ee ee ee eee ~- —- 


ee ete ee er ee re 





TTT PM SR ce yg NE SSG SS 


——— NT FS a o-— 


+ ee ae ee 


C----— 


— ss C—====READ FINITE ELEMENT MESH COORDINATES AND ELEMENT THERMAL 
C——-—---P ARAMETERS 
lee 
ee" READIL, 7150) (XCTPsYC(I), T=l¢_¥M) 
READ(1L,7200) { (NODO( I,J) eJ=1,3) esTKCI)D ,TKT( I) ,TKF(I),ROWC( I), 
LROWCF( I) ,ROWCT( 1) pPCWTR(T) pUNWAT(I),IT=LsNEMS) 
_ (BG EO95 l= isNeis 
 FNOD(I) = 0. 
1005 CMNOD{(T) = 0. 
~~" 1006 TF LBC) 4600+4000,1007 
1007 NSERT = NSURF + 1 
IF ( NSERT .GT. LBC) GO TO L113 
To  READ(19¢7250,END =3800) (NBC(S),BC(J) sJENSERT,LBOC) 
1113 IF { NSURF) LOC8,1CO3,1111 


_ oe -——- — — <i - -— —— oa. -_ — = oe 


-— et eee ee -—— <> ee ee ee >) ee 


ee ee er ee eee we a le es wae — ° ~ 


. 
mw ee ee ee FO ee ees Se ie ee ee ee ores re ——————— eae 


° —-c seem 2 Se = - oe —e ar 
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{ “se @ ee © © eacom aw ow 


= 
—_ ai <i ee en nn -—- — — —_—s - = SOT a Me ET ww | eS — emns conte 


— —  -1111 DOO 1112 T=1lsNSURF 
ll1l2 NBC{I) = I 


(----- 


. 2 ee UNKNOWN NODE POINTS 





mp mm nr er 


1008 IF (TIME) 3900,1009,1355 
1009 I = 1 
~ NKNOD)=1 
1010 DO 1050 J=1,L8C . 
LL=NBC(J) 
IF ( LL-NKN(I) NE. 0) GO TO 1050 > 
NKN(I) = NKNCI) + 1 
1050 CCNTINUE 
—— 1060 1=14) 00 ee SS eee 
IF ( I .GT. AK) GO TO 1070 
NKNCIT)=NKNQT~-1) + 1 


ee ee ee ee ™-e —- Se ee oo - re em ee ees eee 


rs ee 





I ee ee eee 





a ll cA iy crits PMS Sm ee a Se iS re ee 


~ GO TO 1010. 
1070 CONTINUE 
Ue — e 
— C=====COMPUTE THE NUMBER OF UNSPECIFIED NODES COMMCN TO EACH ELEMENT = 
C-~--- 
IF { NSURF .LE. 0 ) GO TO 1300 
— 90.:1225 T=1l,NSURFO 7 ae 
1225 CMNOD(!) = lL. 
NSERT=NSURF 41 
— GOiMTO..1325 ee 


1300 NSERT=1 
1325 DO 135C J=NSERTsNEMS 
~~~ 90 1350 KelyB 0 0 ey SS a a 
L= 0 
DO1L350 M=1,LBC . 
IF (NBC(M)—-NOD(J,K) «EQ. 0) GO TO 1350 Se 
IF (€ NOD(J:K) .€Q. 0 ) GO TO 1350 
L=t+tl : 
— ‘TF ( U .EY. tee) GO TO 13500 id 7 7 i 
CMNOD(J) = CMNOD(J) + 1. 
1350 CONTINUE 


-— ~ 1355 CONTINUE "0 WUT se ae 
IF (NSURF) 1300,1800,1360 
1360 READ(1,7275,END=3800) DATE, TAMB,VEL,CLOUDs,ALTsRHsSRFLXsEVFLX _. 
— ~ [TF (VEL .GT.e Le) GO TO 1380 ae 
VEMP= 1. 


C----- 


oT Cm === ASSIGN APPROPRIATE THERMAL AND ROUGHNESS PARAMETERS | vo ae 
Co---— 
1380 IF ( SNOW) 1385,1385,1395 
— 1385 EMISS = .90 —— ee 
- ALBDG = .15 
GO TO 1430 ; Ae ; 
. 1395 EMISS = .$8 = ae Rar a 
IFC TAMB -— 32.0) 1400,1405,1405 
1400 ALBDGO = .85 
— 60 10 1420 (er ree 
1405 IF (TOAY) 1410,1410,141L5 





~ "3445 IF € SNOW) 1446,1446,1447 


— 2 


_ 1437 CNVCT = 25.#VEL*.38 | 


~“1900 CONTINUE — 
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1410 TOAY = TIME a 
1415 AGE = ABSUTIME - TOAY) 


ALBDO = .62*EXP(-.0215*AGE) + .23*EXP(-.494AGE) 
1420 CONTINUE 


_—— a i a : — Se er — ee ee ee 





ec me mm ee 


C----- 
C—----COMPUTE CONVECTIVE HEAT TRANSFER COEFFICIENT. 

a 

1430 IF ¢ ABSUTAVE "- TAMB) wll. Lei GO TO 40> °° (ooo 


IF (TAVE — TAMB) 1435514401445 
1435 IF (SNOW) 143621436,1437 
1436 CNVCT= 49.6*VEL*.38 
GO TO 1800 





-— a a AS NM AR, SD A EON A Ft ie 


GO TO 1800 
1440 IF ( SNOW) 1441,1441,1442 
1441 CNVCT = 86.6*VEL®.30 

GO TO 1800 a 
1442 CNVCT = 60.*VEL*.30 

GO TO 1800 


LT A 


1446 CNVCT = 120.*VEL¥.22 


IF (NSURF) 1803,1803,1801 
1801 SUM = O. 
~  §$0 1802 I = LeNSURF a eat rit Ul iden 
1802 SUM = SUM + TO (I) 
TAVE = SUM/FLOATUNSURF) 
IF ( ABS(TAVE ~- TAMB) .LT. 100.)} GO TO 1803 
WRITE (3,8200) 
GO To 3900 
1803 IF (TIME) 1805+1805,2020 
1805 NOUM=l1 


{ eee cee EP ae 
rp 


C---—-— 
NCOL=1 BD cps, - 
IF ( KOOE .NE. 129 GO TO 1900 
IFCKODE -EQ. 1) NCOL=2 

GO TO 2010 


~~-oe cm me ee ee ee ee ee ee 


Cm me en em ee rer 


IF(KODE -EQ. 2) GO TO 1950 
NOUM= NASURF + 1 


be a RR I aS 


~ 1950 DO 20C0 N=NDUMeNEMS 


DO 2000 [=l,3 

OO 20CO J=le3 ee ee 
NN = NCO(NeI) — NODCIN,JS) + 1 
2000 IF ¢( NCOL —- NN LT. O 9 NCOL=NN 


= ee — se. eee ee 








—_——— a ce ee 


ee 


te - -_ 2e-=— « . — i 6 $e - . Se a 


eter = = =e 


GO TO 1800 
"1447 CNVCT = 97.*VEL*.22 
GO TO 1800 - 
1800 REDO = 0. 
_ IF(OPT .LT. 1.) GO TO 1810 a 
REDO = le 
1810 CONTINUE 
~ —~NROW = MM an: 
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_.2010 CONTINUE . _ = 
IF ( NCOL .LE. 20 ) GO TO 2020 
WRITE (3,7300) 


____GO TO 3900 
<== nee 
C-—----ZERO SYSTEM ARRAYS 
C--——— 
2020 DO 2040 I=1,MM ao. 0 a 
R(I} = 0. 
__00 2040 J=1,NCOL ees 
~ SIlisd) =e 
S(I,J) = 0 
PtlI,J) = 0. oe 
2040 WIIeJ) = O- a re ee 
IF (NSURF) 20454 2045,2041 
c= 
C—--~COMPUTE INCIDENT SOLAR FLUX AND LONG WAVE ENERGY EXCHANGE 
C—-—-— 


__ 2041 QSOL=(1.-ALBDO)*( 1.-( .82-.024#ALT)*CLOUD/10.)*SRELX _ 
TABS” = (TAMB—-32.)*(5./9.) + 273216 
VALUE = 373.16/TABS 

~SVPL = ~7.902SB¥(VALUE ~1.)# 5.02808*ALOGLO( VALUE) 

1-€1.e3816E-7) =( 10. *¥(11.344%(1.-€1L./VALUE)) )-1.) 

24(8.1328E-3)*(100%*(-3.45149%(VALUE -1e}) ~1.) 

“3 ¢ 3.005715 ae 
SWVP= 10.**(SVPL ) 
VP=SWVP*RH/100. 

“TRANS= 1.0 — .O24*ALT °° °»©— 

ae _ ATMEM = 274 + 29125*CLOUD + .0049 *VvP 
QLRAD={1.—TRANS*CLOUD/10.) #14.389E-7) *(EMESS#( TS#*4. )-ATMEM* 
“LI TABS##4.)) = — 


C-—~--- 


C-—----CCNVERT RADIATION FLUXES TO AN EQUIVALENT TEMPERATURE _ 


eae 


TEQ = (QSOL -— EVFLX - QLRAD) /(.336*CNVCT) 
EQTEM = TAMB + TEQ 
marron pO 2043 {=lsNsuee° «=~ = =+273S © ©6ClhCChCCT eee le 
BCCI) = EQTEM 
2043 TK(I) = .336*CNVCT 
C~—--=- 
C—---CONSTRUCT THE SYSTEM MATRICES 
Cc--- 
2045 DO 2600 M=1,NEMS 
CALL CENTER (M,XSAR,YSAR eNSURF, KODE,DELTY,AM) 
IF €( KODE .EQ. 1) GO TO 2550 
IF (M .LE. NSURF) GU TO 2550 ——— re 
CALL XIETA(M,XBAR, YEAR, XLOC,YLOC) 
CALL AREA (XLOC,YLOC,AM) 
— °° °° °° °° & +f £ &4 .GT. 0.) GO TO 2550 
WRITE(3,57350) M 
GO TO 3900 
SIT oy (OLOI NAT LS al aaa Rs ES 
2050 IF ( M eLEs NSURF ? GO TO 2575 





me eee 





—~- —e  Semeeeee 


—_—_— “cae iF fF == =~ 


Ee TR ee © te 


—— 





EE ee ee See ay ot cee = oo — : eo s+e6 . wee we tees = _— me we ee ee ee ee ee 





—— ee ee ee ee ee ee ee ee ee ee ee ee ee eee ee See ee ee ee - _ 





_ TF ( TIME .GT. 0.) GO TO 2575 
DO 220C K 2 123 
J = 0 
DO 22CC L = 1,t8C 
NCOM = NOD(M,K) 
IF ( NBC(L) — NCOM .NE. 0) GO TO 2075 
GO TO 2200 
2075 IF ( NCCM .EQ. 0 ) GO TO 2200 
J=zJsel 
IF (J .EQ. LBC) GO TO 2100 
GOT TO 2200 ~~ eae ee 
2100 IF ( KCDE .EQ. 1 ) GO TO 2110 
' BETA = AM 
~~~ GO TG 2120, 
2110 BETA = ABS (DELTY)*#AM 
2120 DO 2125 I = 1yNK 
~~ "TE -( NKNCIDY — NCGM .EQ. 0 3 GO TO 2130 
- 2125 CONTINUE 
2130 FLUID(M) = BETA®PCHTROM) SUNWAT(M) *(.20L) 
—< VOL(I) = VOLCI) + SETA/CMNOD(M) _ 


Se el OO tet Oe TT ce ge — ee ew nee eee ——— _ 


A PSF ee 


—— i ee ee eee ee — — ee ee ee epee eee ee ee ee 


ee ee ee ee 


i we we oe) ee oe eee ee -—— —ns- —= - —— ome pa ee eee ee el 





me I eee ee ee 


mm mm mmm em mmr a mm cr er te ee 


WATE(I) = WATECT) # FLUILDOM) /CMNOD(M) 

. RWCT(I) = RWCTCI) + CROWCT(M) *6ETA)/CMNOD(M) 

— s RWCFCT) = RWCECT) + CROWCF(M) *BETA)/CMNOD(M) ~ _ 
If( TOCI) ~GE. TFAZ) GO TO 2200 
FNOD(M) = FNOD(M) + le 


ee mm cm tr re eee 


~ 2200 CGNTINUE 
2575 CONTINUE 
CALL SMALS(MeSEL sAMsAeNSURF,KCDE 2 DELTY) 
CALL SMALP (M,AMePELsNSURFeKODE,DELTY) ~ 
CALL ADDIN (SEL,PEL»M,KODE) 
2600 CONTINUE : 
~~" TF (TIME) 390052610¢2E50_— 
2610 CONTINUE a 
DO 2625 I=12NK 
——" FZHETOCI) = WATEQCI)#HEAT 
RWCFECI) = RwWCFOUI)/VOLCI) 
2625 RWCTCI) = RWCTCI)D/VOLCI) 
em cme es CONTINUE “Seu 06©6[0)06[6[UCGUlUlUlUt™C™C ee 
C——--— 
C——--—-- INSERT SPECIFIED BOUNDARY TEMPERATURES AND MODIFY THE SePs AND 
— Cm——==——R MATRICES ACCCRDINGLY = = OE 


a 


C¢----- 
CALL BNDOIN(NROW,NCOL 2 LBC) 
a nm ee le 
C——----BEGIN TIME DOMAIN SOLUTION (CRANK NICOLSON) 
¢---—- 


Pe SY ee 





~~ —~""90°2700 T=1,NROW- 7 
DO 2700 J=1,NCOL 
WEI, J)= SOCI,J) © (€2-/0T) PCIe J) 
2700 PCI,J) = (2-/0T) *®0P(IT2,J)) — SCI2JS) 
C¢-—~—--- 


C———--PRETRIANGULARIZE THE W MATRIX 


ew 6 ee 


CALL PRESCLONROW, NCOL ) 


ae a 0 em ee ~~ 





mee ee ee 


i  —— ——— eee — re ie oes «© = - 
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—_-=—s—s 2800 CGNTINUE ee 
CALL CCMB(TOs,NROW,NCOL) 
DO 2900 L=lyNROW 
— 2900 G(t) = 2.#R( 1) + Z(T) 
CALL FINSOL (GyNROW,NCOL) 
TIME = TIME + OT 
COUNTI= CCUNTI + l 
ITSEG = LTSEG"¢+ L 
I=1 
oe ee CGO 0.3084 
3081 G(I) = TFAZ 
3082 TO (ft) = Gtl) 
I=[+41 
“3084 CONTINUE © 
IF ( I .GT. NROW) GO TO 36889 
DIFF = G(I) - TFAZ 
IF ( DIFF) 3085,3C081,3487 
C—--- 
C——~TTSIMULATE FREEZING PROCESS 
C—---~ 
3085 IF ( NFREZ(I) -EQ. 1 ) GO TO 3082 
DELQ'= RWCT(I)*VOL(I)*0I1FF 
“FZHET(I) = FZHET(I) + DELQ 
IF ( FZHET(I) .GT. 0.) GO TO 3081 
NFREZ(I) = 1 | 
— GCI) = TEAZ * FZHET(I) ¢% ( RWCFECID*VOLCI) YOU 
WRITE(3,7400) NKN(I) TIME 
FZHET (1) = WATEC(T)*HEAT 


— ee —— — — =~ ~ a ee ee 


um Oe ee Oe ee, 


ee eee wee ae mm Am mm mg ep cep 








te sete ey me Sg I ee ee ee ee ee 


mp km mm mm ee 


_— ———————————e me eS SS ee ee 


GO TO 3400 
¢-—--— 
C——----SIMULATE THAWING PROCESS 
ee ee ae 


3487 IF ( NFREZC(I) NE. 1) GO TO 3C€82 
3500 DELTQ = RWCF(L)*VOL( LT) FOTFFe(-1.0) 


—— i = oe a 


-FZHET( I) = FZHET(I) + DELTQ 
{F (FZHET(I) .GT. O. } GO TO 3081 
NFREZ(I) = 2 
GCI) = TFAZ ~ FZHET(I) 7 ( RWCTCLDEVOL CT 
WRITE(3,7450) AKN(IT),TIME 
FZHET(L) = WATE( I) *HEAT z _ 
3400 NCON = NKN(I) a a a 
060 3415 J=1l,NEMS 
DO 3415 K=1,3 
NCHOS= NOD(J5K) 
IF (« NCHOS ~£Q. NOON) GO TO 3410 
GO To 3415 
— 34@L0.1F ( J.LE.o NSURF) GO TO 3415° 
IF ( NFREZC(I) .£Q. 2 ) GO TO 3412 
FNOO(J) = FNOD(J) + Ll. 
GO TO 3413 oe 
3412 FNOND(J) = FNOD(JS) - le 
C-_—--- 
~ C—=—=-MODIFY ELEMENT THERMAL PARAMETERS WHEN NODAL PHASE CHANGE OCCURS ~~” 


C-_—-— 


EE i A ey a ee mm ey rm ee ee 


me ee ee er ne ee ee 


me rr ee ee eo ee ee et ———— een Ce — oe 8 ee ee Oe a es ee ee oe 


' 
4 
a i ee ee ES = - is mee -—— 8 ee eee ee ee ee eee 


See ee ee Fe ee ee oe eee eee ee ee oc eee oe ef « -—~ @e 
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ee. 3413 TK(J) = ( TKE(S)*FNOD(S) + ( CMNCD(J)-FNOD(J) )*TKTCJ) 2/CMNOD(J) 
ROWC( JS) = ( ROWCF(J)*FNOD(J) + € CMNOD(JS)-FNODIJ) )#ROWCTOS) 2/ 
1CMNOD(J) : 
.___ 3415 CONTINUE 
REOO = 1.0 eo) ee a ae 


GO TO 3082 - 
3689 IF (NSURF) 3694,3654,3690 © 


~¢—-=-= 1 


C——----COMPUTE AVERAGE SURFACE TEMPERATURE 
C—---= 
3690 S00" s oa 
OG 3691 I[=L»NSURF 
__ 3691 SUM = SUM + TO(T) 
' TAVE = SUM/FLOAT(NSURF) — 
QCONV = CNVCT*.336*#(TAVE ~- TAMB) 

ff IF (SNOW) 3694,3654,3692 

3692 IF(TAVE — 32.) 3694,3693,3693...—~—~ 

3693 SNOW = O. 

WRITE(3-8150) TIME 
WR (307525) 000 ce ae eee 
WRITE(3,7550) ( NKNOI)», TOCI) eI=LyNK) © 
WRITE(3*7520) 
~~ 6G TA 3900 
3694 IF (TIME .GE. END) GO TO 3860 
IF (COUNTI «LT. COUNT) GO TO 3695 : 
———  @=—)——i“‘<‘“‘< OC INTE = OW a lL. ne 
WRITE ( 347500) TIME 
WRITE(327520) 
“WRITE(3;71510)". © © °° «0. SS So <5 ie 
WRITE(327520) 
WRITE(3,7550) ( NKNCI)sTOCT)» IT=1eNK) 
WRITE(327520) ~— ‘i se 
IF (NSURF .LT. L) GO TO 3695 
WRITE(3,8050) TIME 
"WRITE(3,8175) CATE. °° °& © © ©&«373»©€©~™COPrrT.. 
WRITE(3,8225) EQTEM,TAMB 
WRITE(3,7520) 
WRITE(3,8000) ~— 
WRITE(347520) 
WRITE(3¢8100) SWVPsQSOL,EVFLX»QLRADsQCONV 
“WRITE(327530) 
3695 IFC(ITSEG «LT. NTSEG) GO TO 3700 
ITSEG = 0 
~~~ “3700 IF (ITSEG) 3710,3710,3720 ~~ 
3710 GO TO 1006 
3720 IF (RECO .GE. 1.) GO TO 1380 _— 
— "~~ GO TG 2800 ~~” 

3800 KRITE (37600) . 
WRITE(3,7500) TIME | oe 
IF (NSURF.LT. 1) GO TO 3801 oS 
WRITE(3,8175) DATE 
3801 WRETE(3s7520) | 

WRITE( 397510) 
WRITE(327520) 


SS een een eee ES 
Pf 


. 
ee ee i ee a a rm ae ee 


+ 4 Seem vee sf cmmoem = 





eee eee 


EN en renee reece erence nnn a 
* 
. 


a A A a | NR SEL ARNG ANG ge 


ee a ee em ee ce ee = Se Ge ee 2 Oe ee 
. 








ee 


_—  — — on ae oe se ee ee ee ee 


eS ee ee apm mmm mem QE 





— ——eeee ee ee ee ee —————_ s == ee eee 








—~— seme we nO eee = 


* 
ee ee ee 





aaa WRITE ( 
WRITE ( 


MNTR= O 
CO 3810 [=lsNK 


0 


381550 
337520) 


3804 MNTR = MNTR + 1 


IF (MNTR eNEo-1) GO TO 38C7 


Se nenten ania a -— = 


WRITE( 38240) 
wRITE( 397525) 
WRITE (378250) 
WRITE(347520) 
3807 CONTINUE 

CK= FZHET(I)/(W 


PRCNT = (1.2 - C 
GO TO 3809 


(_ NKNCT)D, TOCT) » I=L,yNK) meena e — om en eee 


IF (TOCI) - TFAZ) 3810,3804,3810 ~ 


ATE(L)*HEAT) 


KI*100~. 


3808 PRCNT = CK*LOO. 
3809 WRITE(33;8275) NKN(I)»sPRCNT,VOL( I) 


_ 3810 CONTINUE 


) GO To 381l 
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Sa) a ces eee ve ee et ee - . SE an RE FAAP re ge ee 
. 


rm rm i a er ee 


IF ( NFREZ(I) -EQ. 1) GO TO 3808 


oe I a a a a YS Sethe ena 


a IF (MNTR -GTe O 
WRITE(3,8325) 
> 3811 WRITE(3.7530) _ 
3900 CONTINUE = aa 
GO TQ 5000 
- £000 WRITE( 328300) ree cts AS a 
5000 CONTINUE = a 
STOP 
END vee actos ora ache = cs 
| 
=. men ee a ee ee et 
ee = So. +. ee 
{ 


a a A 2s St a a i 2 st ae ee Ge a 
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__ SUBROUTINE CENTER(M,XBAR »YBARyNSURF »KODEsDELTY, AM) 

C THIS SUBROUTINE COMPUTES THE CENTER COORDINATE OF TRIANGULAR 

c ELEMENT M. 
__... COMMON /BLK5/_NOD(15093) __ 

| COMMON /BLK6/ -X(100) 
COMMON /BLK7/ Y(100) 


I = 


— 


J = 


IF (KODE 


IF { 
1K = 


ee ee 


NOD(M,1) _ 
NOD(M,2) _ 
IF ( KODE.EQO. 1) GO TO 2 


M 


NOD(M,3) 


XBAR = ( X(1) + XJ) 


A te a RESP |p A pe A NS 


eEQ. 2 ) GO TO 1 
eLE. NSURF) GO TO 2 


St gat 


+ X(K) ) / 3-6 


+ Y(K) ) /3- 





—)_/ FLOATINSURE) 


YBAR = ( YC(I) + Y(J) 
GO TO 9 
2 DELTY = Y(J) - Y(T) 
AM = l. a 
IF ( KODE «NE. 3 ) GO TO 9° 
AM = ( X(NSURF) - X(1) 
9 CONTINUE 
RETURN 


——END 





SR cr A ee a | mg ee te ete 


—— re 


-_ a PS 


a ee 


ed 


a ee a eee 


oi et ee 


s—— ee ee - ee ee eee ee eee oe 


Ee Se ee 
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| lS pe mg ee cy hg en es Shs rs lp l  SSS> c 
. 


SUBROUTINE XIETA (MyXBAR,YBAR »XLOCyYLOC) 


C____ THIS SUBROUTINE COMPUTES THE LOCAL COORDINATES OF 
C NODES ON TRIANGULAR ELEMENTS M. 
COMMON /BLK5/ NOD(150, 3) = 


_. COMMON /BLK6/ X(100). 
~ COMMON /BLK7/ Y(100) 

DIMENSION XLOC(3),YLOC( 3) 
DO 100 L=l,3 
T= NOD(M,L) 
XLOC({L)= X(1)—XBAR 

100 YLOC(L)= YC(T)— YBAR_ 
RETURN = - 
END 


FP A a ee 





eS RE eh GS A a AE A | ES ES EO EE EE a EE I I oo e.g ees a> clita 
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im SUBROUTINE AREA (XLOC+YLOC,AM) 
C 
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THIS SUBROUTINE CALCULATES THE AREA OF TRIANGULAR 
C ELEMENT Me 
DIMENSION XLOC(3)+YLOC(3)~ 
AM= (XLOC(2)*YLOC(3)—XLOC(3)#YLOC(2)—XLOC(1) ¥YLOC(3) 
— 1 +XLOC(3)*YLOC(1) +XLOC(1)*YLOC(2)-XLOC(2)*YLOC(1) J/20- 
RETURN 
ENO 


—— ee 


a SE PS SS rr A — mo. 
7 


pg A oT) SASS RST Sieh > hy n,m ga gS 


ee ee ee ee —— ~ _— —- ree — ee ee ee ee ee ee Se ae 
mm em mm a 
ee ee a pe am ee 


eR em Ne i GSE ee 


-——ee — 98 ee eee oe ee ee ee ee ee le e+ == ae em 8 ee lee ~~ — 78 — eww © ae em ws SO i ee aS -_ - > ——2o oe eee tee ew 0 ee ee -——« 
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C FOR 


—— SS 
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— — - a . _— - - ~’*. eee ee ee - -_— ——— = cm rr me ee eee ee eee 


SUBROUTINE SMALS(MySEL,AMyA,NSURFYKONE,DELTY)_ 
THIS SUBROUTINE COMPUTES THE SMALL S MATRIX 
TRIANGULAR ELEMENT M, 
COMMON /BLKS/ NOD(150,3) | 
COMMON /BLK6/ X(100) 

COMMON /BLK7/ Y(100) 

COMMON /BLK8/ TK(150) | 
DIMENSION SEL (393)2A(3y 3) 

I = NOD(M)1) 

J = NOD(M,2) ae 

IF ( KODE «EO. 1 ) GO TO 3 : 
IF ( KODE .E€O. 2 ) GO TO 1 

IF_( M «LE. NSURF) GO TO 3 __ 

K = NOD(M,3) 

Al lel) =X( JS) *®Y(K) =X (K) ®Y (9) 
Al(lse2)=X(K )*Y CTP -X CTD *Y(K) 


ee menatiinemenananditie tine ae eae 
© nn pe SSS SS SP a i SS St 


— ee 


ame ee 


A(1s3)=X(1)*¥ (J) -X (J) VY (T) 
Al2s1)= Y(S)-Y(K) 
A(2s2)= Y(K)— Y(T)_ 


Al2s3) = Y(I) =~ YJ) 

A(391) = X(K) -— X(J) 

A(392) = X(T) ~- X(K) pe idaeme eS  8 ee 
A(353) = X(J) ~ X(T) 

DO 2 J = 193 

DO 2 K = 123 


2 SEL(JsK) = Al2,J)#AC2]sK)F(TKI(IM)/(AM#S.) ») 
1 + Al(39J)*A(39K)*( 1 TK(M)/(AM*4.) ) 


api 


a ee eet ee 


7 eae 


5 


a 


GO TO 5 

ELMS = ( TK(M)*AM) /DELTY }}}}©}»}»™©©.OOO 
SEL(1,1) = ELMS#l. 

SEL(122) = ELMS*(~1.) 

SEL (291) = ELMS#(-1.) 

SEL(292) = ELMS#*l. 

00 4 [= 1y3 

SEL(I,3) = 0. ee 6 
CONTINUE 

RETURN 


—— —— a 


END 


ee mm me re em em a 


-_> em ee I ee ee ee ee ee ee 
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oe ~~ _——<_ao —_——-— << ee ee 


t 
____ SUBROUTINE SMALP(M,AM,PEL,NSURF ,KODE,DELTY ) 
C THIS SUBROUTINE CONSTRUCTS A SYMMETRICAL P ELEMENT 
e C SUBSYSTEM MATRIX FOR THE MTH TRIANGULAR ELEMENT 


; _COMMON /BLK9/ ROWC(150) 
DIMENSION PEL(3,3) 
IF ( KODE .EQ. 1 ) GO TO 2 
IF ( KODE .E0. 2 ) GO TO1_ 
IF ( M .LE. NSURF) GO TO 2 
1 EL = ( ROWC(M)*AM) / 12. 
PEL(1l,1) = 2-*EL 
PEL (152)=1.*EL 
PEL (1,3)=1.*EL 
PEL (2,1 )=1 o*EL 
PEL (252)=2 e*EL 
PEL(2,3)=1.*EL 
PEL(3,1)=1.*EL | 
PEL (392)=1.-*EL 
PEL (343)=2 .*EL 


rs ee ee ee 


NE eee 


mmm ea ae mmm nw ee 


SE ee 


0) OS rr 
2 ELP = ( AM*ROWC(M)*DELTY) /6- 
; PEL(1,1) = ELP#2. 
oe CC REE (les) = ELP¥1l. 
PEL(2e1) = ELP#l. 
PEL (2,2) = ELP*2. 
3BONS 1=l1,3 gj 
3 PEL(I,3) = O. 
5 CONTINUE 
RETURN 


. 
cee ee ee 


END 


mmm Sm tt ma th a NS VS 


eee 


eet ee oe ee ea ee 
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A DEMENSTON SEL (3,3) 
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ee ey 


SURROUTINE ADDIN (SEL,PEL,M,KODE) 
COMMON /BLK2/ $( 100,20) 

COMMON /BLK5/ NONN( 150, 3) 

COMMON /BLK3/ P(100,20) __ 


DIMENSION PEL (3,3) 
_IF_( KODE .EQ. 1) GO TO 2. 
IF ( KODE .EQ. 2 ) GO TO1 
IF ( M .LE. NSURF) GO TO 2 

1 00 100 J = 153. 
NR = NOD(MgJ) 
OO 100 K = 1,3 
NC = NOD(M,K) — NOD(MyJ) + 1 


— a ete a et cere I - 


JF -( NC .tE. 0°) GO TO 100 


———"100° CONTINUE — 


S(NR»NC) = S(NR»NC) + SEL (JyK) 
P(NR»NC) = P(NRyNC) + PEL(JsK) 





GO TO 5 
2 DO 200 J=1 2 


— — 


“NR = NOD(MyJ) 

DO 200 K = 1+2 

NC = NOD(MsK) - NOD(Myy) + 1 = 
IF ( NC .LE. 0) GO TO 200 
S( NRyNC) = SO(NR»NC) + SELG UK) 
P(NR»NC) = P(NR»NC) + PEL (JoK) 

2OGecONmnINUE 

5 CONTINUE 


cm ee ue Sm re ee ee 
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eee ee I IS 8 tay SN ie 
= 


cowie em ee ee 


el a oe ee 


eee eee eee ee 
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COMMON /BLK1/ R(100) 
COMMON /BLK2/ S$(100, 
COMMON /BLK3/ P(100, 


é ———— ee 


SUBROUTINE BNDIN (NROW,NCOL,LBC) 
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20) 
20) 


mE SS SS 


COMMON /BLK12/ NBC(50) 
COMMON /BLK13/ B8C(50) 


"C= SET BOUNDARY CONDITIONS IN MATRICES 


C=- 

_0O0 3230 I=1,LBC 
NN=NBC(1) 
NC=NCOL 

3200_NR= NN=NC+1 | 

IF(NR .LT. 1) GO TO 

R(NR)= RO(NR) - SO(NR 
3210 NC=NC-1 

IF(NC .GT. 1) GO TO 

NC=1 

NR=NN 


3220 NR=NR+1 


ae 





“3210 
»NC)*BC(I) 


3200 — 


IF (NR «GT. NROW) GO TO 3230 


NC=NC+1 


(at lt al cS SSS 


IF (NC «GT. NCOL) GO TO 3230 


R(NR)= RINRJ-S(NN,NC)*BC(I) 


___GO TO 3220 
3230 CONTINUE 
C- 


___C=REFORM MATRICES AND ELIMINATE EQUATIONS 


~~ C=AT BOUNDARY CONDITION NODES 


C- 





—_———— 


NN=NBC (1) -I1+1 
NROW=NROW=-1 


IF (NN .GE. NROW) GO TO 3250 


~ DO 3240 NR=NN,NROW 
R(NR)=R(NR +1) 
DO 3240 NC=1,NCOL 
~ S(NR»NC)=S(NR+1,NC) 
3240 P(NRyNC)=P(NR4+1_NC) 
___ 3250 _NR=NN _ 
“N=2 
3260 NR=NR-1 


ae ee 


DO 3300 1=1,LBC _ 


legheett. 1) GO 10 3300:meeee 


N=N+ ] 
IF( N .GT. NCOL) GO 
DO 3270 J=NyNCOL 
P(NRygJ—L)=P (NR od) | 
3270 S(NR»J—LI=S(NRe J) 
P(NR,NCOL )=0.- 
~~ S$ (NR »NCOL)=0. 
GO TO 3260 
3280 S(NR»NCOL)=0. 
— ~ P(NR,»NCOL )=0. 
3300 CONTINUE 
____RETURN_ 
END 


——$S—pamn _——— 





Sa a ee ee Se 


eee ee 


TO 3280 


mS 
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—_—e eee ee 
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_ SUBROUTINE PRESOL (NROW,NCOL) __ 


- C-THIS SUBROUTINE PRE-TRIANGULARIZES A SYMMETRIC 
- C-MATRIC IN BAND FORM FOR SOLUTION BY THE GAUSSIAN 
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—— ee rm mmm wr em erm wr se ee ee ee 


C-ELIMINATION METHOD. FINAL SOLUTION IS BY SUBROUTINE FINSOL. 





COMMON /BLK4/ ST(100,20) 
_ «COMMON. /BLK10/ W(10020) 
N=0 
100 N=N+1 
_ STON, NC OL) SS W (N ? 1) 


> —_. ——— 


= ee 
C-REOUCE PIVOT EQUATION ; 
oe ee 8 
TF(N-NROW) 150,300+ 150 
150 DO 200 K=2»sNCOL 
___ST(N¢K=L)=W(N,K) 
200 W(INsK)= W(NyK)/W(Ny 1) 





C= 

___C-REDUCE REMAINING EQUATIONS 

DO 260 L=2,NCOL 

I=N+t-1 ae. a 
—~ IF(NROW-1I) 26092407240 

240 J=0 
__00 250 K=LeNCOL 
J=J+1 


250 WIJ) SWUT pJ)—-ST(N,L=1) #W(NGK) 
260 CONTINUE 








SSS TI SPS 


i ES 


“G60 TO 100 -_ 
300 RETURN 
_ END fio 


A Go er ER ~ A 


a a oes 
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SUBROUTINE COMB (Y,NROW,NCOL). ce 


= es aE : 
ww C------------—— ------------~-- --- -- -- -- +--+ + + -- -- +--+ +--+ +--+ -- 
peced C- THIS SUBROUTINE MULTIPLIES SYMMETRIC MATRIX P__ 
| C-TO MATRIX Y AND STORES THE RESULT IN 2. oe a - 
G Sam BE Se eee ewe ns 2 SS Sees See eSB ewes eK Sew e es Oe Se wwe ew SWE Oe See wwe ea aeaeawa eee 
C- 





COMMON /BLK3/. P(100,20) 
COMMON /BLK11/ 2Z(100) 
_DIMENSION Y(100) 
DO 200 I=l1,NROW 
Z(IV=YCI)*P(IT9 1) 
DO 200 K = 2,NCOL _ 
L=I-K+l 
IF(L.~LT.1) GO TO 100 
Z(IV=Z(1) + Y(L)*P (LK) 
100 N=I+K-1 -_ 
IF(N.GT.NROW) GO TO 200 
_._ZCY)=Z(1) +Y (NDP (TSK) 
200 CONTINUE 
RETURN 


IR ARR ES 





ee ee eee 
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. ss. = - - - _ —_- _ . or =— 
_—_ = —_— ~~ —< — - _—-— -— oe es oo — -—— a oe —— «a+ » _— ee ee eee ee er ee ee ee ee — —— ee see ee oe ees ee eee ee 


SUBROUTINE FINSOL(SS»sNROW,NCOL) 


— ee —_— = 
oo A ee LS 


C a a a wn wn wn wm wr wn wm wn we or eo eo oe rw we oe oo ee we ee 


et C-THIS SUBROUTINE SOLVES A SET OF LINEAR SIMULTANEOUS 
_ C-EQUATIONS WHOSE COEFFICIENT MATRIX HAS BEFN PRE-TRIANGULARIZED 


| C=BY THE GAUSSIAN ELIMINATIION METHOD. THE SYSTEM MATRIX 
} C=IS IN BAND FORM AND IS SYMMETRICAL. SOLUTION IS PLACED 
_C-IN THE LOAD MATRIX. 7 


f [= Sy Se a ce es cs ce cS es a cee eee ed ee Ge es ee ee eee: —— a a a 
° 
| c= 


— ——— -_ te me os - Se Re re ae 


___COMMON /BLK4/ $T(100,20) 
COMMON /BLK10/ W(100,20) 
DIMENSION SS(100) 


——— — oe 


~ €=REDUCE LOAD VARIABLES 
C- 

N=0 

went!  _— ——_  \o . 
SS(N) = SS(N)/ST(NyNCOL) 
IF (N-NROW) 110,200,110 

~ 210 CONTINUE a a 

DO 130 K=2,NCOL 
=N+K=-1 
~~ JFFC(NROW-L) 130,120,120 
120 SS(L)I=SS(LIAKST(NyK-1)#SS(N) 


—— wee 





C= 

_C-BACK SUBSTITUTION _ 

C- 

200 N=NROW 

300 N=N-1 

~ - “JFIN) 350,500,350 

350 DO 400 K=2,yNCOL 
L=N+K-1 : 
TF(NROW-L) 40053707370 ©”. ” © ©) 

370 SS(N)= SS(N) — WIN SK)*SS(L) 

: 400 CONTINUE 

— sco 90 300 " °°» 
500 RETURN 

END 


* 
se ies we eee em a ee ee ee eee 
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PS pe 


N379Yd AVWILYVd Juv S3IGON ONIMONIOd 341 


re | re eee eee ee =~ I 


A ee a A ee i ee ee ee i ee ee ee SS = ee a eee ee eS ee ee me A ee Se eee 


em ei eee et ee ——e ee oe a ———t oe ee Hee we se we eee ae ee 


O¥°ZE “@2 O¥°ze L2 Zo°2E 92 »°2e G2 £#8b*°2€  °&42-b2  °&2€S$°2€  °©€2 
QS°ZE 22 S9°ZE TZ TL°2E 02 B2°ZE 61 G8°ZE 81 26°2E | L6°ZE 9T 
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